^ in suspended bilayer graphene: the interplay of disorder and band gap 



O 
(N 

> 

O 



CD 



a 



o 

> 

O 
O 



o 



X 



D. S. L. Abergcl, H. MinQ E. H. Hwang, and S. Das Sarma 

Condensed Matter Theory Center, Department of Physics, 

University of Maryland, College Park, MD, 20742, USA 

We present an interpretation of recent experimental measurements of ^ in suspended bilayer 
graphene samples. We demonstrate that the data may be quantitatively described by assuming a 
spatially varying inter-layer potential asymmetry (which generates a band gap) induced by local 
electric fields resulting from charged impurity disorder in the graphene environment. We demon- 
strate that the fluctuations in the inter-layer potential asymmetry and density vary between different 
samples, and that the potential asymmetry fluctuations increase in magnitude as the density fluc- 
tuations increase. This indicates that the mechanism causing this effect is likely to be extrinsic. 
We also provide predictions for the optical conductivity and mobility of suspended bilayer graphene 
samples with small band gaps in the presence of disorder. 



Recently, the fabrication of suspended bilayer graphene 
(BLG) devicesiS, has allowed the direct local measure- 
ment of the thermodynamic quantity X = -j^ using 
scanning single electron transistor spectroscopy. This 
is an important quantity to study since it is linked di- 
rectly to thermodynamic observables such as the ther- 
modynamic density of states, electronic compressiblility 
and the quantum capacitance of the interacting electron 
liquid in the BLG. In samples mounted on a Si02 sub- 
strate, the level of disorder is generally too high for the 
intrinsic physics near the charge neutrality point to be 
seen because of the existence of strong inhomogeneity 
(electron- hole puddles) induced by the disorder"^. How- 
ever, suspended samples provide a setting where the ef- 
fective strength of the impurities is reduced as reflected, 
for example, in the enhanced mobility^i^ and it is likely 
that the low-density physics is more readily accessible. 
In the experiments^'-^, a region of decreased compress- 
ibility (or equivalently a peak in K) is seen near the 
charge-neutrality point in zero magnetic field indicating 
the possible presence of a small gap in the low energy 
band structure. It is claimedi that the size of the gap 
which generates this peak is approximately 2nieV, and 
that the magnitude of charge inhomogeneity is of the or- 
der of 10^°cm~^. The authors of Ref. |l| state that Castro 
et alM- show that this level of in-plane charge inhomogene- 
ity by itself would create a gap which is a factor of 10 
smaller than that observed, ruling out disorder-induced 
fluctuations as a mechanism for the generation of this 
gap. Instead, they claim^ that this peak is evidence for a 
correlated ground state of the interacting electron liquid, 
such as one of those proposed in the recent literature^r— . 
Some of these proposed states may induce a spontaneous 
many-body charge transfer instability in the bilayer sys- 
tem which could open a small band gap at zero density, 
and the region of decreased compressibility seen in the 
experiments may be evidence for the formation of one of 
these charge-transfer-instability statesi. However, we be- 
lieve that this interpretation of the experimental results 
is problematic because the out-of-plane charge imbalance 
which would drive a disorder-induced gap is not the same 
quantity as the inplane charge inhomogeneity which Mar- 



tin et al. usei to estimate the gap size. Therefore it is 
important to look more closely at the role of disorder in 
the suspended bilayer graphene system to investigate in 
depth the issue of the possible existence (or absence) of 
a many-body quantum phase transition creating a spon- 
taneous band gap. The goal of the current work is to 
carry out such a phenomenological investigation to see 
whether random charge impurities in the bilayer envi- 
ronment could produce a signature in the compressibility 
consistent with the observations of Refs. [l| andy. 

Transport measurements suggest that there are some 
residual impurities in the environment of the suspended 
graphene, each of which will have an electric field associ- 
ated with it, and it is well established that the electrons in 
the BLG order themselves into puddles of electrons and 
holes in order to screen these fields^. Additionally, it 
was shown recently ^2, that these electric fields also cause 
a spatial variation in the interlayer potential asymmetry 
which is highly correlated with the disorder profile. This 
is a different effect from the inplane reorganization of 
charge, and the two effects may in principle be of quite 
different magnitudes, and both are likely to be present 
in disordered bilayer graphene. Although this interlayer 
effect was reported for bilayer graphene supported on a 
substrate^^, the general principle applies to suspended 
samples as well. 

In this article, we revisit the issue of density fluctu- 
ations induced by charged impurity disorder, and show 
that our calculation of K by averaging over the density of 
the electrons in the puddles^ can be applied in two dif- 
ferent phenomenological ways to this situation and gives 
results for the gap fluctuations, density fluctuations, and 
required charge imbalance which are self-consistent. By 
fitting this theory to measurements of K from two dif- 
ferent devicesi^, we show that the gap and charge inho- 
mogeneity are device-dependent, which is, of course, not 
unexpected since disorder in different devices is likely to 
be different leading to different interlayer potential asym- 
metry in each sample. We also demonstrate that the 
transport characteristics associated with our phenomeno- 
logical disorder fits are reasonable, and predict that in an 
optical spectroscopy experiment, the gap will still be ob- 



scured by disorder in these samplesi^. Since our results 
are internally consistent, we claim that the mechanism 
of local density fluctuations induced by charged impurity 
disorder cannot be ruled out as a cause of the observedi'^ 
peak in K at low carrier density. Our work by no means 
eliminates the possible occurence of a many-body BLG 
instability, it just establishes the qualitative and quanti- 
tative importance of disorder in thinking about the quan- 
titative aspects of the observed effects, and points out 
that there is a possible single-particle extrinsic cause for 
the observed effective gap, namely, the inter-layer poten- 
tial asymmetry induced by the random charged impuri- 
ties in the BLG environment. In reality, both interac- 
tion and disorder may be present in a non-perturbative 
manner, making a microscopic theoretical analysis quite 
challenging, particularly since the graphene landscape 
becomes inhomogeneous due to the formation of electron- 
hole puddles in the presence of the charged impurity dis- 
order. 

In order to model the clean BLG system, we consider 
the continuum limit of the tight-binding Hamiltonian. 
In this approach, the wave functions are described by 
Bloch functions modulated by a spinor, the elements of 
which correspond to the four sites in the BLG unit cell^. 
We make this assumption despite the fact that the inho- 
mogeneity in the disorder potential breaks translational 
symmetry because we assume that the changes in the 
potential landscape are spatially slowly varying enough 
that the physics associated with p-n junctions or quan- 
tum dots do not manifest. We also point out that this 
theory has been very successful in describing capacitance- 
based compressibility experiments in samples which are 
much dirtier than those we currently discussiS. The band 
structure is determined by the relative strength of elec- 
tron hopping between the four lattice sites. To simulate 
the experimental situation as accurately as possible, we 
go beyond the standard nearest-neighbor tight binding 
formalism (where vp is the Fermi velocity of monolayer 
graphene which comes from the intra-layer hopping be- 
tween adjacent atoms and 71 is the energy associated 
with the inter-layer dimer bond which generates the fi- 
nite effective mass at low energy) by including the next- 
nearest neighbor inter-layer hops with energy 74 and the 
onsite energies given by the lattice site asymmetry A and 
a potential asymmetry u between the two layers. This 
potential asymmetry directly generates a band gap at 
the K point of the band structure, and we shall there- 
fore use the terms 'inter-layer potential asymmetry' and 
'gap' synonymously in this paper. The Hamiltonian for 
one valley is then 
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where U4 = v374a/2?i with a the lattice constant, and 
k and are respectively the radial and angular part of 



the wave vector measured from the K point. The spec- 
trum is valley-symmetric and isotropic. The effect of 
the 74 and A terms is to introduce an asymmetry be- 
tween the low-energy conduction and valence bands, and 
has previously been used to explain an observed asym- 
metry in capacitance measurementsSi. Using the results 
collated in Ref. [20 we choose representative parameters 
71 = 0.4eV, 74 = 0.15eV, A = O.OlSeV, vp = lO^ms"! 
and a ~ 0.246nm. We describe the band structure and 
K in the absence of disorder in the appendix. 

The existence of puddles of carriers in graphene is 
a well-established phenomenon both in experiment and 
in theory22. In a previous paper^^, we introduced the 
macroscopic averaging over density fluctuations and ap- 
plied this theory to measurements of the capacitance of 
BLG devices finding good agreement with the experimen- 
tal data. The idea of the theory is that the potential fluc- 
tuations associated with charged impurity disorder gener- 
ates a spatially- varying on-site term in the tight-binding 
Hamiltonian. This in turn leads to a local fluctuation in 
the Fermi energy and hence in the carrier density. The 
statistics of this fluctuation can be encapsulated within a 
distribution function P. The coupling between the tip of 
the SET and the electron liquid in the bilayer graphene 
is capacitative, so we assume that each region with a 
given carrier density has a local capacitance associated 
with it. Then, the total capacative coupling of the SET 
to the sample is the parallel coupling of all of the areas, 
which equates to the geometrical average over the region 
sampled by the SET. Therefore, the areal average can be 
replaced by an average over the charge density distribu- 
tion P. We note that the dimension of the region sampled 
by the SET in the experiments of Ref. [l| is ~ lOOnm 
whereas the puddle size seen in Ref. [ij is ~ lOnm so 
that the averaging should be reasonable in this case. We 
stress that this philosophy deliberately retains the inho- 
mogeneity of the system and treats it as the primary 
manifestation of disorder. This is a fundamentally differ- 
ent approach from, for example, a diagrammatic expan- 
sion of the electron-impurity interaction which contains 
within it an average over realizations of disorder which 
explicitely restores the translational symmetry. We be- 
lieve that the important disorder physics to be included 
in the low-density regime of graphene near the charge 
neutrality point is the formation of electron-hole pud- 
dles in the system around the discrete quenched charged 
impurity centers leading to strong spatial inhomogene- 
ity modelled by the distribution function P mentioned 
above which has been calculated in the literature in a 
few instances^ and measured in experiment for a bare 
Si02 substrate^^. 

For density fluctuations associated with puddles of car- 
riers in graphene systems we write the local density as 
n{r) = nc + h{r) where Uc is the average density induced 
by external gates and the fluctuations are encapsulated 
within h. We parameterize these fluctuations by their 
standard deviation Sn, and computing the average of K 
over the probability distribution for the density results 



in the following expression for the average inverse com- 
pressibility 
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The distribution function P is talken as the following 
Gaussian 
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since there is compelling theoretical- and 
experimental^^ evidence that this is approximately 
the correct form. In this theory, the layer asymmetry 
u — Uc is constant in space (hence the subscript 'c') and 
is generated by either external electric fields (such as 
those from gates) or by an asymmetry of charge between 
the two layers of the BLG generated by the intrinsic 
electron-electron interactions. Analytical evaluation of 
this procedure is not possible since, for finite 74 and A, 
the eigenvalues of the Hamiltonian ([1]) are the roots of 
a quartic polynomial. Therefore, the expression for /i 
as a function of n, its derivative, and the integration in 
Eq. ([2]) must be calculated numerically. 

Figure [ija) shows the results of finding the best fit be- 
tween Eq. ^ and the experimental dataii^. We find 
that the published data^ (which we shall denote by the 
label PRL) is quantitatively described by a constant gap 
Uc = 1.7meV with density fluctuations characterized by 
6n = 0.8 X 10^*'cm^2, which is very close to the val- 
ues claimed by Martin et al. in their paperJ^. The un- 
published data^ (denoted by PC, as in 'private com- 
munication') is described by Uc = 3.3meV and 6n = 
2.1 X 10^*'cm^2. Therefore, a different degree of in-plane 
charge homogeneity and different size of gap is required 
to describe the experimental data in each sample which is 
perfectly reasonable since we expect the disorder to have 
sample-to-sample variations. Also, the band gap required 
to fit the data is larger for the more strongly disordered 
sample. The opposite would be true if the gap was the 
result of an intrinsic interaction-driven effect since in this 
case, the gap due to interactions would be universal, but 
disorder would work to reduce its sizei^. Therefore, these 
results, particularly the fact that the more strongly dis- 
ordered sample exhibits the larger effective experimental 
band gap, indicate that there is some role being played 
by extrinsic effects in these experiments. 

We now introduce another step in this theory, moti- 
vated by recent observations of spatial fluctuations in 
the interlayer potential asymmetry in bilayer graphene^^. 
The interlayer asymmetry was measured in electron pud- 
dles and in hole puddles of a bilayer graphene flake which 
was mounted on an Si02 substrate. The difference in the 
asymmetry of approximately 70meV was seen between 
the two puddles at zero back gate voltage, and an as- 
sociated band gap was seen in the spectrum of Landau 
levels measured via STM. Although we expect that the 
magnitude of this effect will be smaller in suspended sam- 
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FIG. 1. (a) Fits of the density averaged K from Eq. ([2]) to the 
experimental data. Note that this theory includes a spatially 
uniform band gap Uc. The black lines are for published data^, 
and the red lines for the unpublished data^ . (b) Fits of the 
double averaged K from Eq. (Q with Uc = to the experi- 
mental data. Note that in this theory, the band gap u = u{r) 
fluctuates in space, and the distribution of these fluctuations 
is parameterized by 5u. The table shows a comparison of the 
fitting parameters for the two samples and two procedures. 



pies, the presence of disorder- induced puddles (as demon- 
strated above by our fitting to the appropriate theory) 
means that it is sensible to assume that the same disor- 
der may also induce (albeit smaller in magnitude) local 
fluctuations in the potential asymmetry just as it does in 
graphene on substrates in Ref. 17. In this case, the poten- 
tial asymmetry becomes a function of position such that 
u = u(r) = Uc + u{r) . The fluctuations are contained 
within u, and we assume that they can be described by 
the distribution P with standard deviation Su (just as 
the density fluctuations can) because the same disorder 
is generating both types of fluctuation. Note that it is the 
in-plane component of the electric fleld associated with 
disorder which causes the formation of puddles, while 
the spatially varying inter-layer potential asymmetry is 
driven by the out-of-plane component of the electric field. 
Therefore, while we expect the two types of fiuctuation to 
be correlated, it is not obvious exactly what relationship 
should exist between them. However, on a phenomeno- 
logical level, we can treat dn and the potential asymmetry 
fiuctuations parameterized by Su as independent param- 
eters. We can therefore perform an average over the gap 
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FIG. 2. Plots of distributions P and Q for (a) xq < 5x, (b) 
a;o = 5a;, and (c) xq > 5a;. In all cases, we have 5x = 1. 



fluctuations in the same way we previously did for the 
density. The spatially varying part of the band gap u 
is described by 5u, and the Gaussian is centered around 
Uc- The spatial average then corresponds to the following 
average over the distribution function: 

A'(ne,(5n, Uc, (5u) = K{nc,Sn,u)P{u,Uc,Su)du. (4) 

We call this procedure 'double averaging'. In order to as- 
sess if the observed gap might be entirely due to disorder- 
induced fluctuations, we fit this expression with Uc = to 
the experimental data. Therefore the fitting parameters 
are now Su and Sn, so we have two parameters as was 
the case for the density averaging procedure. The results 
are displayed in Fig. [IJb). We find that the PRL data 
is fitted by Su = 2.2meV and Sn = 0.75 x lO^^cm^^ and 
the PC data by Su = 5.2meV and Sn = 3.2 x lO^Ocm^^ 
The fit of the double-averaged theory is actually slightly 
better than that which is averaged only over density (es- 
pecially in the region of the peak), but both averaging 
procedures give results which are in reasonable agree- 
ment with the experiment and therefore it is not possible 
to distinguish between a spatially fluctuating disorder- 
induced gap and a uniform intrinsic gap in the current 
data^i. Therefore, our work convincingly demonstrates 
they key importance of disorder, particularly the density 
and potential fluctuations associated with the inhonio- 
geneous puddles, in determining the experimental com- 
pressibility measurements of Refs. HI andy. 

Using the parameters from these fits, we can make the 
following observation. In principle, if there was truly no 
disorder then both Sn and Su should be zero. However, if 
we take the two pairs extracted from the double averag- 
ing fitting procedure and assume a linear extrapolation 
to Sn = 0, we find a residual Su ~ l.SmeV. This may in- 
dicate the presence of an intrinsic (possibly, many-body) 
gap not generated by the fiuctuations, so we examine that 
issue now. The experimental measurements of K in the 
PRL data suggest that the size of the uniform gap, if it 
exists, is approximately 2meV. Since this is of the same 
order as the gap predicted by extrapolating to the clean 
limit, we have run fits of our double-averaged theory with 
a constant gap of this size to the experimental data. We 



FIG. 3. Optical conductivity for (a) Uc = 5meV and (b) 
Uc = 2meV for various degrees of disorder broadening. The 
color of the lines is the same in both panels. Dashed lines 
indicate the clean case. 



find that there is no improvement in the fitting for any 
value of Sn or Su. To explain this, we need to analyze the 
distribution function for the gap fiuctuations. We know 
that K is an even function of the gap size u, and there- 
fore the absolute value |u| is the key for determining the 
distribution of K. Therefore, although we average over 
the Gaussian P in Eq. (j4]), the effective distribution of 
K is non-Gaussian, and given by 

Q{x, xq, Sx) = P{x, xo, Sx) + P{—x, xo, Sx), (5) 

where x > and P{x,xq,Sx) is defined in Eq. ^. Fig- 
ure [2] shows the distribution Q for the cases xq < Sx, 
xq = Sx, and xq > Sx. We see that in all but the last 
case, the distribution is dominated by small values of x 
because of the significant contribution from the P{—x) 
term, and it is only when Sx > xq that the Gaussian 
shape is revealed in Q. In the two sets of data which we 
have studied, we see that Uc found from the single av- 
eraging procedure and the Su from the double averaging 
obey the relation Uc < Su. This suggests that the disor- 
der is too strong to determine if a spontaneous gap exists 
in the absence of disorder. In order to measure the spon- 
taneous gap, samples with levels of disorder such that 
Su < Uc are needed. This necessitates the development of 
samples with much higher quality (i.e. less disorder), per- 
haps mobilities which are a factor of five or so larger, to 
decisively establish the existence of a spontaneous many- 
body BLG gap. In short, the disorder-induced gap needs 
to be a factor of 2 or more smaller than the many-body 
gap for an unambiguous conclusion on this issue. 

The gap size may also be extracted from measurements 
of the optical conductivity via absorption or reflection 
experiments^^"— . To demonstrate that this may be an 
unreliable process in the case of the small band gaps (of 
the order of a few meVs) we are discussing in this con- 
text, we plot in Fig. |3]the optical conductivity derived 
from the Kubo formula with u = 5meV in panel (a) 
and with u = 2meV in panel (b). By way of compar- 
ison, in the linear (high density) part of the hyperbolic 
spectrum of bilayer graphene, the density fluctuations Sn 
correspond to fluctuations in the Fermi level of approx- 
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gives an energy scale of approximately 20meV. When 
the Fermi energy is in the sombrero (low density) region 
with u — 5meV, the flatness of the bands means that this 
density fluctuation corresponds to a fluctuation of less 
than ImeV. We include the disorder through a constant 
broadening of the bands parameterized by the energy T 
which appears in the Green's function as 
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The optical conductivity is calculated from this Green's 
function from the following expression^*^: 
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FIG. 4. (a) Calculated mobility of suspended bilayer graphene 
as a function of impurity density for different strengths of the 
short-ranged impurities, (b) Conductivity calculated within 
effective medium theory as a function of gate voltage Vg for 
a charged impurity density, n^ = 2 x lO^^cm"^, and a fixed 
short range potential ndVjf = 0.3 (eVA)^. Difi^erent curves 
correspond to the band gap Eg = 0, 0.15, 0.3, 0.33 eV (from 
top to bottom). Inset shows the calculated conductivities as a 
function of carrier density for two different charged impurity 
densities m = 2 x lO^^cm"^ (dashed line) and rii = 4 x 
lO^^cni"^ (solid line) with zero band gap. 



Figure [3] shows that for F = ImeV, the peak due to 
the onset of the intraband conductivity is already signif- 
icantly blurred. By the time F > Uc, the peak is com- 
pletely smeared, as shown by the F = 5meV line. We 
mention in this context that the same conclusion can also 
be reached by considering an inhomogeneous broadening 
effect rather than a homogeneous broadening as used in 
Eq. ([6|) - for example, we can introduce a density fluc- 
tuation leading to an uncertainty in the Fermi energy by 
ImeV, leading to exactly the same conclusion that the 
optical gap measurement would be unable to discern the 
small band gap in the presence of this inhomogeneous 
broadening effect. 

We also predict the transport properties of suspended 
bilayer graphene with disorder so that independent char- 
acterization of the samples in which K was measured 
can be done and compared with our theory incorporat- 
ing inhomogeneous puddle effects of density fluctuations. 
We use the highly successful Boltzmann-RPA formalism, 
which has been used to study the transport in mono- 
layer graphene and BLG systems^^. This theory incor- 
porates both long-range Coulomb impurities and short- 
range scatterers. In Fig. Hl^a) we show the mobility of 
suspended BLG as a function of charged impurity den- 
sity rii, calculated at a carrier density n = 10-^^cm^^ for 
different short-ranged disorder potentials, UdV^ , where 
Ud is the density of short-ranged impurities and Vq is the 
strength of the potential. If we assume a charged im- 
purity density n^ sa 5n, then we have /.t~4 — 5x10^ 
cm^/Vs with a fixed short range potential naV^ — 0.3 
(eVA)^. However, we note that even though the car- 
rier density fluctuation is related to the impurity density, 
the theoretical relation between two densities is not yet 
known. In Fig. HJb) we show the conductivity calculated 
within effective medium theory— as a function of gate 



voltage Vg in the presence of an energy band gap. We 
assume that the energy dispersion of carriers is given by 
Eki, = vh?k'^ /{2m) + vEg/2, where v ~ ±1 indicates the 
electron [v — \) or hole {v = —1) band, m is the effective 
mass, and Eg is the band gap. In this figure a charged im- 
purity density n.^ = 2 x 10^''cm^^ and a fixed short range 
potential rirfVJ]^ = 0.3 (eVA)^ are used. Our calculation of 
Fig. mb) shows that due to the electron hole puddles the 
transport gap (i.e. the region of zero conductivity) does 
not occur even though there is a band gap (top three 
curves in the figure). Only large enough energy gaps 
(the lowest curve) reflect the transport gap. However, 
our calculation also shows that for small density fluctu- 
ations (i.e. shallow electron-hole puddles) the transport 
gap can be induced by the smaller energy gap. In the 
inset to Fig. S^b) the conductivity of suspended BLG 
calculated within the effective medium theory^ is shown 
for two different charged impurity densities, rii = 2 and 
4xl0^°cm^^ for a flxed short range potential n^VQ = 0.3 
(eVA)^. Our calculation shows the minimum conductiv- 
ity at the charge neutral point is higher for a low mobility 
(more strongly disordered) sample. 

In conclusion, we have demonstrated that even though 
suspended graphene samples are known to be signifl- 
cantly less disordered than those mounted on substrates, 
the effect of disorder in recent measurements of -# can- 
not be ruled out. Using phenomenological averaging 
procedures, which specifically incorporate the effects of 
strong inhomogeneous broadening arising from electron- 
hole puddle induced spatial density and potential fiuctu- 
ations, we have demonstrated that current experimental 
data cannot distinguish between an intrinsic, spatially 
uniform band gap generated by electron-electron inter- 
actions, and a spatially fluctuating band gap induced by 
disorder. Specifically, these disorder-induced fiuctuations 



can account for the observed data without invoking the 
presence of an intrinsic gap^. In order to be confident 
about the origin of the band gap, higher quality samples 
where u^ > Su must be fabricated. One possibility is to 
carry out measurements on BLG samples fabricated on 
BN substrates which typically have much lower charged 
impurity disorde r^^'^^ . 

The theory we develop involves exactly one uncon- 
trolled approximation, assuming the density fluctuations 
and the fluctuations in u to be uncorrclatcd. Going be- 
yond this approximation would require knowledge about 
the details of the local electric fields producing the fluc- 
tuations in n and u which is currently unavailable ex- 
perimentally. It would be straightforward to include the 
correlation in our theory if microscopic information about 
the correlator between fluctuations in n and u is available 
from experiments or some other microscopic calculations. 
In the absence of such information we assumed them to 
be independent. 

We also emphasize that the theory presented in this 
article is phenomenological and is a first step toward the 
understanding of the role of disorder in the experiments 
of Refs.[l| and[3. What we have achieved in this paper is a 
demonstration that spatial density and potential fluctua- 
tions are capable of producing a spatially fluctuating gap 
consistent with the experimental data. Many questions 
of details could be raised with respect to our phenomeno- 
logical procedure, for example, the averaging procedure 
we employ is certainly not unique, but in the presence of 
real spatial fluctuations which explicitly break the trans- 
lational invariance in graphene near the Dirac point, our 
procedure is physically motivated as it averages precisely 
over the quantities which fluctuate according to some 
distribution function which has been calculated in the 
literature2*i^i^>^. A better (purely numerical) method 
to tackle the problem would be to use the techniques 
developed in Refs. y, [la, and |32 to explicitly and self- 
consistently calculate the real-space electronic structure 
in the presence of the electron-hole puddles generated 
by the random charge impurities, and then to use this 
real-space electronic structure to directly numerically cal- 
culate the compressibility. Such a calculation is obvi- 
ously numerically prohibitively difficult, and our work is 
a short-cut (albeit a rather simple one) in carrying out 
such a numerical procedure. Since the Gaussian prob- 
ability distribution function we use is numerically well- 
verified in Refs. y, [la, [22, and [32| through self-consistent 
calculations, we believe that our results, although phe- 
nomenological, should be qualitatively and quantitatively 
valid. The eventual theory involving a microscopic calcu- 
lation of the compressibility in the presence of long-range 
Coulomb disorder and electron-electron interaction in- 
cluding the realistic band structure of bilayer graphene 
is very far in the future since no one knows how to ap- 
proach such a non-perturbative problem even at a formal 
level. 
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FIG. 5. (a,b) The band structure of clean bilayer graphene 
with a band gap. The dashed lines correspond to the electron- 
hole symmetric band structure found with 74 = A = [see 
Eg. lAlj and the solid line to the full band structure with 74 = 
150meV, A = 18meV. Black lines are for u = 20meV and red 
lines for u = 60meV. (c) The bottom of the conduction band 
for u = 20meV showing the absence of sombrero structure 
for finite 74, A. (d) The bottom of the conduction band for 
u — 60meV showing the restoration of the sombrero structure 
for finite 74, A. (e) ^ for clean bilayer graphene. The lines 
correspond to the same as in (a). 
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Appendix A: Clean bilayer graphene 

The band structure associated with the Hamiltonian 
presented in Eq. ([T]) is shown in Fig. [S] In the symmetric 
case with 74 = A = 0, the low energy bands have the 
dispersion 



E±=± 



:<±+n^{jf + u^) (Al) 



where k — hvpk. Writing this as a function of density 
and differentiating yields the following analytical expres- 
sion for -J^ in this limit in the conduction band: 

an 



1-2 2 



K 



2 V^^^+^aA^+^ 

h^l-TT l-(M^+7?)/(2a) 



A < u^, 
A> w^ 



(A2) 



^+T + ^-" 



where X = h VpTrn and a 



A(w2 +72) + :^. In con- 



trast, when 74 and A are included in the analysis, analyt- 
ical solutions of the eigenvalue problem are not possible. 
The energy spectrum is found as the roots of the follow- 



ing quartic polynomial equation for E: 



E^ - 2AE^ 



A^-7?-y-2.^(l + X^ 



Au2 



2 

7? - A^) 



2k2(A + x'A + 2x7i) 

2 



-«Mi-x^) 



^" 



i; 



= (A3) 



where x = v^/vp- The Fermi surface may be ring-shaped 
for finite values of u so that the wave vectors correspond- 
ing to the inner and outer Fermi surfaces k± are given by 
the equation nn = k'^ — k'^ with the constraint that k-^- 
and k- must give the same energy for the appropriate 



band when substituted into Eq. (jA3l) . When the Fermi 
energy is above the sombrero region, fc_ = 0. The ex- 
tracted value of fc+ can be substituted into Eq. (|A3[) to 
find /i. The differentiation to obtain K can the be car- 
ried out numerically. Evaluations of these equations are 
shown in Fig.[5]for the band structure and K for two dif- 
ferent values of the band gap, and with and without the 
tight-binding parameters 74 and A. The step feature in 
K is a, result of the change in topology of the Fermi sur- 
face from a ring to a disc as the Fermi energy leaves the 
sombrero region. For small values of u, the conduction 
band does not show the sombrero shape when 74, A are 
finite and this step is absent. The spike at the origin is a 
(5-function resulting from the discontinuity of the Fermi 
energy at zero density as it jumps from the valence band 
to the conduction band. 
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